{
 "metadata": {
  "signature": "sha256:99a159af5f25d7a29b5a367bf7eb812a7ed122f2fe4a942270743af3f4c35177"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Solving large Markov Chains\n",
      "===========================\n",
      "\n",
      "This page shows how to compute the stationary distribution pi of a large\n",
      "Markov chain. The example is a tandem of two M/M/1 queues. Generally the\n",
      "transition matrix P of the Markov chain is sparse, so that we can either\n",
      "use scipy.sparse or Pysparse. Here we demonstrate how to use both of\n",
      "these tools.\n",
      "\n",
      "Power Method\n",
      "------------\n",
      "\n",
      "In this section we find pi by means of iterative methods called the\n",
      "Power method. More specifically, given a (stochastic) transition matrix\n",
      "P, and an initial vector pi\\_0, compute iteratively pi\\_n = pi\\_{n-1} P\n",
      "until the distance (in some norm) between pi\\_n and pi\\_{n-1} is small\n",
      "enough.\n",
      "\n",
      "Fist we build the generator matrix Q for the related Markov chain. Then\n",
      "we turn Q into a transition matrix P by the method of uniformization,\n",
      "that is, we define P as I - Q/l, where I is the identity matrix (of the\n",
      "same size as Q) and l is the smallest element on the diagonal of Q. Once\n",
      "we have P, we approximate pi (the left eigenvector of P that satisfies\n",
      "pi = pi P) by the iterates pi\\_n = pi\\_0 P\\^n, for some initial vector\n",
      "pi\\_0.\n",
      "\n",
      "More details of the above approach can be found in (more or less) any\n",
      "book on probability and Markov Chains. A fine example, with many nice\n",
      "examples and attention to the numerical solution of Markov chains, is\n",
      "\\`Queueing networks and Markov Chains' by G. Bolch et al., John Wiley,\n",
      "2nd edition, 2006.\n",
      "\n",
      "You can get the source code for this tutorial here:\n",
      "[tandemqueue.py](files/attachments/Solving_Large_Markov_Chains/tandemqueue.py)\n"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#!/usr/bin/env python\n",
      "\n",
      "labda, mu1, mu2 = 1., 1.01, 1.001\n",
      "N1, N2 = 50, 50\n",
      "size = N1*N2\n",
      "\n",
      "from numpy import ones, zeros, empty\n",
      "import scipy.sparse as sp\n",
      "import  pysparse\n",
      "from pylab import matshow, savefig\n",
      "from scipy.linalg import norm\n",
      "import time"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This simple function converts the state (i,j), which represents that the\n",
      "first queue contains i jobs and the second queue j jobs, to a more\n",
      "suitable form to define a transition matrix."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def state(i,j):\n",
      "    return j*N1 + i"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Build the off-diagonal elements of the generator matrix Q."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def fillOffDiagonal(Q):\n",
      "    # labda\n",
      "    for i in range(0,N1-1):\n",
      "        for j in range(0,N2):\n",
      "            Q[(state(i,j),state(i+1,j))]= labda\n",
      "    # mu2\n",
      "    for i in range(0,N1):\n",
      "        for j in range(1,N2):\n",
      "            Q[(state(i,j),state(i,j-1))]= mu2\n",
      "    # mu1\n",
      "    for i in range(1,N1):\n",
      "        for j in range(0,N2-1):\n",
      "            Q[(state(i,j),state(i-1,j+1))]= mu1\n",
      "    print \"ready filling\""
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "In this function we use scipy.sparse"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def computePiMethod1():\n",
      "    e0 = time.time()\n",
      "    Q = sp.dok_matrix((size,size))\n",
      "    fillOffDiagonal(Q)\n",
      "    # Set the diagonal of Q such that the row sums are zero\n",
      "    Q.setdiag( -Q*ones(size) )\n",
      "    # Compute a suitable stochastic matrix by means of uniformization\n",
      "    l = min(Q.values())*1.001  # avoid periodicity, see the book of Bolch et al.\n",
      "    P = sp.speye(size, size) - Q/l\n",
      "    # compute Pi\n",
      "    P =  P.tocsr()\n",
      "    pi = zeros(size);  pi1 = zeros(size)\n",
      "    pi[0] = 1;\n",
      "    n = norm(pi - pi1,1); i = 0;\n",
      "    while n > 1e-3 and i < 1e5:\n",
      "        pi1 = pi*P\n",
      "        pi = pi1*P   # avoid copying pi1 to pi\n",
      "        n = norm(pi - pi1,1); i += 1\n",
      "    print \"Method 1: \", time.time() - e0, i\n",
      "    return pi"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now use Pysparse."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def computePiMethod2():\n",
      "    e0 = time.time()\n",
      "    Q = pysparse.spmatrix.ll_mat(size,size)\n",
      "    fillOffDiagonal(Q)\n",
      "    # fill diagonal\n",
      "    x =  empty(size)\n",
      "    Q.matvec(ones(size),x)\n",
      "    Q.put(-x)\n",
      "    # uniformize\n",
      "    l = min(Q.values())*1.001\n",
      "    P = pysparse.spmatrix.ll_mat(size,size)\n",
      "    P.put(ones(size))\n",
      "    P.shift(-1./l, Q)\n",
      "    # Compute pi\n",
      "    P = P.to_csr()\n",
      "    pi = zeros(size);  pi1 = zeros(size)\n",
      "    pi[0] = 1;\n",
      "    n = norm(pi - pi1,1); i = 0;\n",
      "    while n > 1e-3 and i < 1e5:\n",
      "        P.matvec_transp(pi,pi1)\n",
      "        P.matvec_transp(pi1,pi)\n",
      "        n = norm(pi - pi1,1); i += 1\n",
      "    print \"Method 2: \", time.time() - e0, i\n",
      "    return pi"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Output the results."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def plotPi(pi):\n",
      "    pi = pi.reshape(N2,N1)\n",
      "    matshow(pi)\n",
      "    savefig(\"pi.eps\")\n",
      "if __name__ == \"__main__\":\n",
      "    pi = computePiMethod1()\n",
      "    pi = computePiMethod2()\n",
      "    plotPi(pi)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Here is the result:\n",
      "\n",
      "![](files/attachments/Solving_Large_Markov_Chains/pi.png)\n",
      "\n",
      "## Improvements of this Tutorial\n",
      "\n",
      "Include other methods such as Jacobi's method or Gauss Seidel. "
     ]
    }
   ],
   "metadata": {}
  }
 ]
}